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Abstract 

In  this  paper,  we  develop  nonlinear  constitutive  equations  and  resulting  system  models  quan¬ 
tifying  the  nonlinear  and  hysteretic  field-displacement  relations  inherent  to  lead  zirconate  titanate 
(PZT)  devices  employed  in  atomic  force  microscope  stage  mechanisms.  We  focus  specifically  on 
PZT  rods  utilizing  (I33  motion  and  PZT  shells  driven  in  d^i  regimes,  but  the  modeling  framework  is 
sufficiently  general  to  accommodate  a  variety  of  drive  geometries.  In  the  first  step  of  the  model  de¬ 
velopment,  lattice-level  energy  relations  are  combined  with  stochastic  homogenization  techniques  to 
construct  nonlinear  constitutive  relations  which  accommodate  the  hysteresis  inherent  to  ferroelectric 
compounds.  Secondly,  these  constitutive  relations  are  employed  in  classical  rod  and  shell  relations  to 
construct  system  models  appropriate  for  presently  employed  nanopositioner  designs.  The  capability 
of  the  models  to  quantify  the  frequency-dependent  hysteresis  inherent  to  the  PZT  stages  is  illustrated 
through  comparison  with  experimental  data. 
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1  Introduction 


Stage  mechanisms  employing  the  ferroelectric  material  lead  zirconate  titanate  (PZT)  have  played  a 
fundamental  role  in  scanning  tunneling  microscope  (STM)  and  atomic  force  microscope  (AFM)  design 
since  their  inception  due  to  the  high  set  point  accuracy,  large  dynamic  range,  and  relatively  small 
temperature  sensitivity  exhibited  by  the  compounds  [10].  To  illustrate,  consider  the  prototypical 
AFM  design  depicted  in  Figure  1.  To  ascertain  the  3-D  surface  structure  of  a  sample,  it  is  moved 
laterally  along  a  pre-determined  x-y  grid  by  a  PZT-driven  stage.  The  response  of  a  highly  flexible 
microcantilever  to  changing  atomic  surface  forces  is  monitored  by  a  reflected  laser  beam  measured 
via  a  photodiode,  and  forces  corresponding  to  the  cantilever  displacement  changes  are  determined 
via  Hooke’s  law.  A  feedback  law  is  used  to  determine  voltages  to  a  transverse  PZT  stage  which 
produces  displacements  in  the  z-direction  to  maintain  constant  forces.  A  complete  scan  in  this 
manner  provides  a  surface  image  of  the  compounds.  Additionally,  PZT  actuators  are  often  used 
to  drive  the  microcantilevers  at  resonance  to  achieve  the  tapping  mode  operation  used  to  reduce 
damage  to  specimens. 

Two  representative  stage  designs  are  depicted  in  Figure  2.  The  first  employs  stacked  PZT  actua¬ 
tors  utilizing  (I33  electromechanical  motion  to  achieve  longitudinal  positioning  along  the  pre-specified 
x-y  grid.  A  second  stage  provides  the  transverse  motion  required  to  ascertain  the  sample  topography. 
Rod  models  with  linear  and  nonlinear  electromechanical  input  relations  are  constructed  to  quantify 
the  PZT  transducer  dynamics  in  this  design.  The  second  geometry  employs  a  cylindrical  shell 
with  half  poled  ^33  to  provide  horizontal  (x-y)  motion  and  half  poled  ^31  for  vertical  (z)  motion 
as  depicted  in  Figure  2(b)  —  to  enhance  vibration  isolation  and  reduce  hysteresis  and  constitutive 
nonlinearities.  Thin  shell  models  are  developed  to  characterize  this  stage  design. 

To  illustrate  issues  which  must  be  addressed  by  models,  field-displacement  data  from  the  stacked 
actuator  depicted  in  Figure  2(a)  is  plotted  in  Figures  3  and  4.  The  data  in  Figure  3  was  collected  at 
0.1  Hz  and  illustrates  the  nested,  hysteretic  relation  between  input  fields  and  generated  displacements 
in  a  nearly  quasistatic  regime.  The  data  in  Figure  4  was  collected  at  frequencies  ranging  from  0.279  Hz 
to  27.9  Hz  to  illustrate  the  frequency-dependence  of  the  hysteresis  as  well  as  certain  dynamic  effects. 

At  low  frequencies,  the  inherent  hysteresis  can  be  accommodated  through  proportional-integral- 
derivative  (PID)  or  robust  control  designs  [5,  6,  18,  24].  However,  at  the  higher  frequencies  required 
for  applications  including  real-time  monitoring  of  biological  processes  (e.g.,  protein  unfolding),  com¬ 
prehensive  product  diagnostics,  and  single  electron  spin  detection  [23,  34],  increasing  noise-to-data 


Figure  1:  (a)  Configuration  of  a  prototypical  AFM,  and  (b)  surface  image  determined  by  one  lateral 
sweep. 
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Figure  2:  Actuator  configurations  employed  for  sample  positioning  in  AFM:  (a)  stacked  actuators 
employed  as  x-  and  y-stages,  and  (b)  cylindrical  PZT  transducer. 


ratios  and  diminishing  high-pass  characteristics  of  control  filters  preclude  a  sole  reliance  on  feedback 
laws  to  eliminate  hysteresis. 

Alternatively,  it  is  illustrated  in  [16,  17],  that  use  of  charge-  or  current-controlled  amplifiers  can 
essentially  eliminate  hysteresis.  However,  this  mode  of  operation  can  be  prohibitively  expensive  when 
compared  with  the  more  commonly  employed  voltage-controlled  amplifiers,  and  current  control  is 
ineffective  if  maintaining  DC  offsets  as  is  the  case  when  the  re-stage  of  an  AFM  is  held  in  a  fixed 
position  while  a  sweep  is  performed  with  the  y-stage. 

The  need  to  significantly  increase  scanning  speeds  with  general  amplifiers  motivates  the  devel¬ 
opment  of  models  and  model-based  control  designs  which  accommodate  the  frequency-dependent 
hysteresis  inherent  to  the  PZT  actuators  employed  in  the  AFM  stages.  As  detailed  in  [25],  there 
exist  a  number  of  general  approaches  and  frameworks  for  quantifying  the  constitutive  nonlinearities 
and  hysteresis  in  the  general  class  of  ferroelectric  materials  which  encompass  PZT.  These  include 
phenomenological  macroscopic  models  [20],  Preisach  models  [9,  22],  domain  wall  models  [28,  29], 
nricronrechanical  models  [4,  14,  15],  mesoscopic  energy  relations  [3,  13]  and  homogenized  energy 
models  [27,  33].  Within  the  context  of  AFM  design,  Croft,  Shed  and  Devasia  [5]  have  employed  a 
combination  of  a  viscoelastic  creep  model  and  nonlinear  Preisach  representation  to  compensate  for 
hysteresis  and  creep  in  an  AFM  stage  whereas  a  domain  wall  model  was  employed  in  [30]  for  the 
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Figure  3:  Nested  minor  loops  in  0.1  Hz  field-displacement  data  from  a  stacked  PZT  stage  of  the  type 
depicted  Figure  2(a). 
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(a)  (b)  (c) 

Figure  4:  Frequency-dependent  field-displacement  behavior  of  a  stacked  PZT  stage  of  the  type 
depicted  in  Figure  2(a):  sample  rates  of  (a)  0.279  Hz,  (b)  5.58  Hz,  and  (c)  27.9  Hz. 


characterization  of  hysteresis  in  certain  stage  constructs.  Primary  requirements  for  nonlinear  hys¬ 
teresis  models  for  the  PZT  actuators  in  an  AFM  are  (i)  flexibility  with  regard  to  frequency-dependent 
hysteresis  effects  —  the  frameworks  of  [5,  30]  are  limited  in  this  regard  — ,  (ii)  exact  or  approximate 
invertibility  for  linear  control  design,  and  (iii)  sufficient  efficiency  for  real-time  implementation  at 
the  speeds  required  for  present  and  future  applications. 

In  this  paper,  we  develop  AFM  transducer  models,  based  on  a  homogenized  energy  framework  for 
characterizing  hysteresis  and  constitutive  nonlinearities  in  ferroelectric  materials,  which  meet  these 
criteria.  On  Section  2,  we  summarize  the  framework  developed  in  [33]  for  quantifying  hysteresis 
in  the  field-polarization  relation  and  develop  constitutive  equations  which  characterize  the  elastic 
and  electromechanical  behavior  of  the  PZT  material.  These  constitutive  relations  are  employed  in 
Section  3  to  construct  rod  and  shell  models  for  the  stages  depicted  in  Figure  2,  and  the  well-posedness 
of  the  models  is  established  in  Section  4.  Numerical  approximation  techniques  for  the  transducer 
models  are  summarized  in  Section  5,  and  the  capability  of  the  framework  to  quantify  the  biased 
and  frequency-dependent  hysteresis  behavior  of  the  transducers  is  illustrated  in  Section  6  through  a 
comparison  with  the  experimental  data  plotted  in  Figures  2  and  3. 

With  regard  to  criteria  (ii)  and  (iii),  the  construction  and  experimental  implementation  of  model 
inverses  to  linearize  the  nonlinear  dynamics  is  demonstrated  in  [12].  Hence  the  models  provide  an 
effective  framework  for  characterizing  the  hysteresis  and  nonlinear  dynamics  inherent  to  PZT-based 
nanopositioners  in  a  manner  which  promotes  stage  and  control  design. 

2  Constitutive  Relations 

In  this  section,  we  summarize  the  development  of  constitutive  relations  which  quantify  the  nonlinear 
and  hysteretic  map  between  input  fields  E  and  stresses  a  and  the  polarization  P  and  strains  e  gen¬ 
erated  in  ferroelectric  materials.  These  relations  are  developed  in  three  steps.  In  the  first,  Helmholtz 
and  Gibbs  energy  relations  are  constructed  at  the  lattice  level  to  quantify  the  local  dependence  of 
P  and  e  on  E  and  a  for  regimes  in  which  relaxation  due  to  thermal  processes  is  either  negligible  or 
significant.  In  the  second  step  of  the  development,  material  nonhomogeneities,  polycrystallinity,  and 
variable  field  effects  are  incorporated  through  the  assumption  that  certain  material  properties  are 
manifestations  of  underlying  distributions  rather  than  constants.  Stochastic  homogenization  in  this 
manner  yields  macroscopic  models  which  quantify  the  bulk  hysteretic  E-P  behavior  measured  in  fer¬ 
roelectric  materials.  Finally,  necessary  conditions  associated  with  minimization  of  the  Gibbs  energy 
are  invoked  to  obtain  1-D  and  2-D  constitutive  relations  quantifying  the  elastic  and  electromechanical 
behavior  of  the  transducer  materials. 


3 


2.1  Helmholtz  and  Gibbs  Energy  Relations 

As  detailed  in  [33],  an  appropriate  Helmholtz  energy  relation  is 

e)  =  i/jp(P)  +  ^Ye2  -  aieP  -  a2eP 2  (1) 

where  the  component 

'  \n{p  +  PR)2  ,  p<-Pi 

i/jp(P)  =  <  hn(P  -  Pr )2  ,  p>  Pi 

_  (Pr  -  PR)  (g  -  Pr)  ,  |P|  <  Pi 

quantifies  the  internal  energy  due  to  dipole  processes.  As  shown  in  Figure  5,  Pj  is  the  positive 
inflection  point  which  delineates  the  transition  between  stable  and  unstable  regions,  Po  denotes  the 

unstable  equilibrium,  and  Pr  is  the  value  of  P  at  which  the  positive  local  minimum  of  V’  occurs. 

The  parameter  r/  is  the  reciprocal  of  the  slope  of  the  E-P  relation  after  switching  occurs.  The 
second  term  on  the  right  side  of  (1)  quantifies  the  elastic  energy  whereas  the  third  and  fourth  terms 
quantify  electromechanical  coupling  effects.  Here  Y  denotes  the  Young’s  modulus  and  ai,a2  are 
electromechanical  coupling  coefficients. 

The  Gibbs  energy  relation 

G(E,  a,  P ,  e)  =  ipP{P)  +  ^Ys2  -  ai eP  -  a2eP 2  -  EP  -  ae  (2) 

incorporates  the  elastic  work  ae  and  electrostatic  work  EP.  This  provides  the  functional  that  is 
minimized  or  balanced  with  the  relative  thermal  energy  to  provide  local  E-P  relations  and  global 
electromechanical  constitutive  equations.  The  reader  is  referred  to  [25,  27]  for  details  regarding  the 
manner  through  which  the  Gibbs  energy  incorporates  the  dependent  variables  e  and  P  in  terms  of 
the  independent  variables  a  and  E. 


V(P)=G(0,P) 


(b) 


Figure  5:  (a)  Helmholtz  energy  if)  and  Gibbs  energy  G  for  a  =  0  and  increasing  fields  E.  (b)  Switch 
in  the  local  polarization  P  that  occurs  as  E  is  increased  beyond  the  local  coercive  field  Ec  given  by 
(6)  in  the  absence  of  thermal  activation. 
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2.2  Polarization  Kernel  —  Negligible  Thermal  Activation 


For  operating  regimes  in  which  thermal  activation  is  negligible,  the  local  E-P  relation  is  determined 
from  the  necessary  conditions 

dG  n  d2G  n 
dP  ~  ’  HP2  >  ' 

For  the  piecewise  quadratic  functional  (2),  this  yields  a  polarization  kernel  of  the  form 


P(E)  = 


E 


+  5 


PrV + Saie 


(3) 


r )  —  2a,2£  r/  —  2ci2£ 

where  5  =  1  for  positively  oriented  dipoles  and  5  =  —  1  for  those  having  negative  orientation.  To 
specify  5,  and  hence  P,  more  specifically  in  terms  of  the  initial  dipole  orientations  and  previous 
switches,  we  employ  Preisach  notation  and  take 


[P(E-Ec,0](t)  =  { 


[P(E-,Ec,Om  ,  r(t)  =  0 

E  D~  ,  r(t)  /  0  and  F(maxr(t))  =  —  Ec 

,  r(t)  /  0  and  F(maxr(t))  =  Ec. 


—  —  Pr 

T]  -rt 


(4) 


l  %  +  Pr 


Here 


'  -  -Pr 

T]  ■n' 


[P(f?;f;c,O](0)= 


f  +  Pr 


,  P(0)  <  —  Ec 
,  -Ec  <  E( 0)  <  £c 
,  P(0)  >  Fc 


(5) 


defines  initial  kernel  values  in  terms  of  the  parameter  ^  =  4^  ±  Pr,  0  designates  the  empty  set,  and 
the  set  of  switching  times  is  given  by 

r(t)  =  {ts  G  (0,  t]  |  E(ts)  =  -Ec  or  E(ts )  =  Ec}. 

The  local  coercive  field 


Ec  =  v{Pr  ~  pi) 


(6) 


quantifies  the  field  at  which  the  negative  well  ceases  to  exist  and  hence  a  dipole  switch  occurs.  To 
illustrate,  the  condition  r  ^  0  and  P(maxr(i))  =  Ec  designates  that  switching  has  occurred  and  the 
last  switch  was  at  Ec;  hence  the  local  polarization  is  [P(E\  Ec,  £)](i)  =  +  Pr. 

Remark  1  For  the  drive  levels  employed  for  nanopositioning,  the  stress  effects  on  the  polarization 
are  typically  negligible  which  motivates  taking  £  =  0  in  (3)-(5).  Hence  the  relations 

1 


P(E)  =  -E  +  Pr5 
V 


[p(E-Ec,om  =  i 


[P(E-,Ec,Om  ,  r(t)  =  0 


(7) 


li-P* 


,  r(t)  /  0  and  F(maxr(t))  =  —  Ec 
,  r(t)  /  0  and  F(maxr(t))  =  Ec 


and 


r  E 


[p{E-Ec,om  =  { 


—  —  Pr 

T)  P1 


{  %  +  pR 

are  usually  employed  when  characterizing  AFM  stages. 


E( 0)  <  —Ec 
—Ec  <  E{ 0)  <  Ec 
E( 0)  >  Ec 
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2.3  Polarization  Kernel  —  Thermal  Activation 

If  thermal  relaxation  is  significant,  the  Gibbs  energy  G  and  relative  thermal  energy  kT /V  are  bal¬ 
anced  through  the  Boltzmann  relation 

p{G)  =  Ce~GV/kT.  (8) 

Here  k  is  Boltzmann’s  constant,  V  denotes  a  reference  volume  chosen  to  ensure  physical  relaxation 
behavior,  and  C  is  chosen  to  ensure  integration  to  unity  for  the  complete  set  of  admissible  inputs. 
As  detailed  in  [25,  33],  this  yields  the  local  polarization  relation 


P  =  x+  ( P+ )  +  X-  {P-) . 


(9) 


The  fractions  x+  and  x_  of  positively  and  negatively  oriented  dipoles  are  quantified  by  the  differential 
equations 

x+  =  -p+_x+  +p_+x_ 


which  can  be  simplified  to 
through  the  identity 


X_  =  —  p (_x_  +  p. I x_|_ 

X+  =  -P+-X+  +  P-+(1  -  X+) 


x_|_  +  x_  =  1. 

The  expected  polarizations  due  to  positively  and  negatively  oriented  dipoles  are 

r-Pi 


/  Pe-G(E,P)v/kTdp 

(P+)=J~fh -  -  (p~)  = 


pe-G(E,P)V/kTdp 


-G(E,P)V/kTdp 


—Pj 


-G{E,P)V/kTdp 


(10) 


'Pi 


where  the  denominator  results  from  the  evaluation  of  C  in  (8).  The  likelihoods  of  switching  from 
positive  to  negative,  and  conversely,  are  given  by 


rPi 


e-G{E,P)V/kTdp 


i — Pi+e 


e-G(E,P)V/kTdp 


P+-  = 


tPi-e 


T(T) 


-G{E,P)V/kTdp 


P-+  = 


I- Pi 


T(T) 


i — Pi+e 


-G(E,P)V/kTdp 


(11) 


>Pi-e 


where  e  is  taken  to  be  a  small  positive  constant.  The  relaxation  time  T  is  the  reciprocal  of  the 
frequency  at  which  dipoles  attempt  to  switch.  It  is  proven  in  [25,  33]  that  P  given  by  (9)  converges 
to  the  local  polarization  (7)  in  the  limit  kT /V  — >  0  of  negligible  thermal  activation. 


Remark  2  In  (10)  and  (11),  we  use  the  notation  G(E,P)  to  indicate  that  we  take  £  =  a  =  0  in  (2) 
in  accordance  with  the  assumption  that  stress  effects  on  the  polarization  are  negligible  at  the  drive 
levels  employed  in  AFM  stages. 
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2.4  Macroscopic  Polarization  Model 


The  local  polarization  relations  (7)  and  (9)  exhibit  the  behavior  depicted  in  Figure  6  and  provide 
reasonable  characterization  of  the  E-P  behavior  of  certain  single  crystal  compounds.  However, 
to  incorporate  the  effects  of  material  and  stress  nonhomogeneities,  poly  crystallinity,  and  variable 
effective  fields  Ee  =  E  +  Ej,  we  assume  that  the  interaction  field  Ej  and  local  coercive  field  Ec 
given  by  (6)  are  manifestations  of  underlying  distributions  rather  than  constants.  If  we  designate 
the  associated  densities  by  zq  and  zq,  the  macroscopic  field-polarization  behavior  is  quantified  by  the 
relation 


[p(E)m 


/: 


zq  {Ec)u2{E!)  [P(E  +  Er,  Ec,  0]  (t)  dET  dEc 


(12) 


where  the  kernel  P  is  given  by  (7)  or  (9). 

As  detailed  in  [25,  27],  the  densities  zq  and  zq  are  assumed  to  satisfy  the  physical  criteria 


(i)  zq  (x)  defined  for  x  >  0, 

(ii)  v2(—x)  =  v2(x), 

(hi)  \MX)\  <  cie~aiX , 

\u2{x)\  <  c2e~a^ 


(13) 


for  positive  ci,  a±,  c2,  a2-  The  restricted  domain  in  (i)  reflects  the  fact  that  the  coercive  field  Ec  is 
positive  whereas  the  symmetry  enforced  in  the  interaction  field  through  (ii)  yields  the  symmetry 
observed  in  low- field  Rayleigh  loops.  Hypothesis  (iii)  incorporates  the  physical  observation  that  the 
coercive  and  interaction  fields  decay  as  a  function  of  distance  and  guarantees  that  integration  against 
the  piecewise  linear  kernel  yields  finite  polarization  values. 

Approximation  of  (12)  through  Gaussian  quadrature  techniques  yields  the  approximate  relation 

Ni  Nj  ~ 

[P(E)\(t)  =  J2T,^ECi)MEIj)\P(EIj  +  E]ECi,Zj)}(t)viwj  (14) 

t=i  j= i 

where  E j. .  Ec%  denote  the  abscissas  associated  with  respective  quadrature  formulae  and  Vi,Wj  are 
the  respective  weights  —  e.g.,  see  [25].  Algorithms  used  to  implement  the  approximate  polarization 
model  (14)  are  provided  in  [25,  33]. 

Techniques  for  identifying  the  densities  zq  and  zq  are  illustrated  in  [25,  27].  For  certain  appli¬ 
cations,  reasonable  accuracy  is  provided  by  a  priori  functions  satisfying  the  physical  criteria  (13) 
and  having  a  small  number  of  parameters  to  be  estimated  through  least  squares  fits  to  data  —  e.g., 
variances  and  means  in  normal  and  lognormal  relations.  For  more  general  applications  requiring  high 
accuracy  for  a  wide  range  of  operating  conditions,  the  Nt  +  Nj  discretized  density  values  zq  (ECz)  and 
^2 (Eij)  can  be  estimated  through  least  squares  techniques. 


Figure  6:  Hysteron  provided  by  (a)  the  relation  (7)  with  negligible  thermal  relaxation,  and  (b)  the 
relation  (9)  which  incorporates  relaxation  mechanisms. 
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Remark  3  From  the  perspective  of  both  numerical  and  experimental  implementation  and  the  estab¬ 
lishment  of  the  well-posedness  of  resulting  transducer  models,  it  is  important  to  quantify  the  regularity 
between  input  fields  and  the  polization  predicted  by  (12).  In  Appendix  A,  it  is  established  that  P  given 
by  (12)  is  continuous  with  respect  to  E. 


2.5  Constitutive  Relations 


To  obtain  a  elastic  constitutive  relations,  the  equilibrium  condition 


is  invoked  to  obtain 

o  =  Ye  —  a\P  —  a2  P2 

which  reduces  to  Hooke’s  law  when  P  =  0.  To  incorporate  internal  damping,  we  posit  that  in  the 
absence  of  electromechanical  effects,  stress  is  proportional  to  a  linear  combination  of  strain  and  strain 
rate  (Kelvin- Voigt  damping  hypothesis).  Finally,  we  note  that  the  PZT  stage  mechanisms  are  poled 
and  hence  operate  about  the  remanence  polarization  P  =  Pr  rather  than  the  depoled  state  P  =  0. 
When  combined  with  the  polarization  model  (12),  this  yields  the  1-D  constitutive  relations 


u  =  Ye  +  Ci  -  ai{P  -  PR)  -  a2(P  -  PR )2 

p  oo  poo 

[P(E)](t)=  /  /  miEME^PiE  +  E^E^OmdE^E, 

J  0  J — oo 


where  C  is  the  Kelvin- Voigt  damping  coefficient.  These  relations  are  employed  when  constructing 
rod  models  to  characterize  the  hysteretic  dynamics  shown  in  Figures  3  and  4  for  the  stacked  actuators 
employed  in  the  stage  construction  depicted  in  Figure  2(a). 

The  constitutive  behavior  of  the  PZT  shell  depicted  in  Figure  2(b)  differs  from  that  of  the  rod  in 
two  fundamental  aspects:  (i)  the  longitudinal  actuation  is  due  to  d^i  rather  than  ^33  mechanisms, 
and  (ii)  longitudinal  and  circumferential  stresses  and  strains  are  coupled  due  to  the  curvature.  To 
designate  the  coupled  material  behavior,  we  let  ex,ox  and  £g,og  respectively  denote  the  normal 
strains  and  stresses  in  the  longitudinal  and  circumferential  directions  and  we  denote  shear  strains 
and  stresses  by  exo  and  oxg.  Finally,  we  let  v  denote  the  Poisson  ratio  for  the  material.  The  resulting 
2-D  constitutive  relations 

Y  c  1 

&x  =  1 - 9  (£x  +  V£g)  +  - - 9  (ex  +  nig)  —  - -  [oi(-P  —  Pr)  +  02 (P  —  Pr )2] 

1  —  v-  1  —  V-  1  —  UL  J 


(16) 


/*oo  p  00 

[P(E)](t)  =  /  /  ^(EME^iPiE  +  EYE^OWdEjdE, 

J  0  J — 00 


Y  C  1 

1 - ~{£9  +  u£x)  +  1 - 7^(£9  "b  u£x)  —  y - -  [al(-P  —  Pr)  +  a2 (P  ~  Pr)2] 


1  —  v2 
Y 


1  —  vl 
C 


Gxe~  2(1  +  nfx9  +  2(1  +  nfxe 


are  employed  when  constructing  transducer  models  for  cylindrical  nanopositioning  stages. 


3  Transducer  Models  for  Stacked  and  Cylindrical  AFM  Stages 


We  now  employ  the  1-D  constitutive  relation  (15)  and  2-D  relation  (16)  to  construct  models  for  the 
stacked  and  cylindrical  AFM  stages  depicted  in  Figure  2.  For  the  stacked  actuator,  we  consider  two 
frameworks:  (i)  a  distributed  PDE  model  which  quantifies  displacements  along  the  rod  length  as  a 
function  of  the  input  field,  and  (ii)  a  lumped  model  which  exploits  the  assumption  of  uniform  stresses 
and  fields  along  the  rod  length  to  motivate  an  ODE  quantifying  displacements  only  at  the  rod  end. 
A  comparison  between  characterization  capabilities  provided  by  the  two  frameworks  is  provided  in 
Section  6.  For  the  cylindrical  shell  design,  we  summarize  a  Donnell-Mushtari  model  which  quantifies 
vertical  motion  provided  by  the  ^-component  of  the  stage  depicted  in  Figure  2(b). 

3.1  Rod  Model  for  the  Stacked  Actuator 


Distributed  Rod  Model 


We  consider  first  the  development  of  a  distributed  rod  model  which  quantifies  the  displacement 
u(t,  x)  along  the  rod  length.  In  accordance  with  present  stage  design,  one  end  of  the  rod  is  assumed 
fixed  while  the  other  encounters  resistance  due  to  the  connecting  mechanisms.  We  assume  that 
this  latter  contribution  can  be  modeled  as  a  damped  elastic  system  with  mass  rri( ,  stiffness  kg  and 
damping  coefficient  eg.  The  density,  cross-sectional  area  and  length  of  the  rod  are  denoted  by  p,  A 
and  £  and,  in  accordance  with  (15),  the  Young’s  modulus  and  Kelvin- Voigt  damping  parameter  are 
denoted  by  Y  and  C. 

Force  balancing  yields  the  relation 


d2u  _  dAf 
^  dt2  dx 


(17) 


where  the  resultant  N  =  fA  adA  is  given  by 


c)n  fP“v 

N  =  YAd. x  +  CAd^dt  ~  ai[P{E)  ~  Pr]  ~  a2[P{E)  ~  Pr]2 

once  the  linear  relation  e  =  ^  is  employed  for  the  strains  in  (15).  The  nonlinear  and  hysteretic  map 
between  input  fields  E  and  the  polarization  P  is  specified  by  (12).  The  fixed-end  condition  yields 
u(t,  0)  =  0  and  balancing  forces  at  x  =  £  yields  the  energy  dissipating  end  condition 


fill  fP‘1! 

N(t,  l)  =  - kLu(t ,  £)- cL  —  (: t ,  i)  -  Ml  (: t ,  i) . 

Finally,  initial  conditions  are  taken  to  be  u(0,x)  =  uq(x)  and  |^(0,x)  =  u\(x).  This  provides  a 
strong  formulation  of  the  stacked  actuator  model. 

To  define  a  weak  or  variational  form  of  the  model  which  is  appropriate  for  well-posedness  analysis, 
approximation,  or  control  design,  states  z  =  (u(-),u(£))  are  considered  in  the  state  space  X  = 


I  ci 

x=0  x= / 


Figure  7:  Rod  of  length  £  and  cross-sectional  area  A  with  a  fixed  end  at  x  =  0  and  energy  dissipating 
boundary  conditions  at  x  =  l. 
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L2( 0,£)  x  M  with  the  inner  product 


(<f>i,  <h2)x  =  /  pAfafodx  +  me<pi<p2 
Jo 


(18) 


where  4>i  =  (<?!>i,  ipi),  <h2  =  (<j> 2,  P 2 )  with  p\  =  (j>i(£),  p 2  =  0i(^).  The  space  of  test  functions  is  taken 
to  be 

V  =  {$>  =  ((j),<p)  e  X  \(j)  e  H\ 0,  £),  0(0)  =  0,  0(£)  =  </>} 
with  the  inner  product 

ft 

($i,$2)y  =  /  Y A(l)\(f)'2dx  +  kcpi<p2-  (19) 

Jo 

Multiplication  by  4>  e  i7g(0,£)  =  {<^>  G  i?1(0,f)  |  0(0)  =  0}  and  integration  by  parts  in  space  yields 
the  weak  model  formulation 


.  d2u 

pA—~y(pdx  + 

Jo  L 


. .  ,  du  d2u 

A&c  +  A  drift 


^  dx=  f  fJpdx 
ax  Jo 


+^4  [ai(P  -  PR)  +  a2(P  -  Pr)2]^  ^dx  - 
which  must  be  satisfied  for  all  cj)  £  V. 


du ,  ,  ,  d2 


(20) 


u . 


fau(t,£)  +  C£—(t,l )  +  me-Qp(t,i 


m 


Lumped  Rod  Model 

The  assumption  that  fields  and  stresses  are  uniform  along  the  rod  length  motivates  the  conclusion 
that  strains  (relative  displacements)  also  exhibit  negligible  x-dependence.  Since  the  position  of  the 
sample  is  dictated  by  the  position  of  the  rod  tip  at  x  =  £,  this  motivates  the  development  of  a  lumped 
model  which  quantifies  ut(t)  =  u(t,£). 

From  the  assumption  of  uniform  strains  along  the  rod  length,  we  take 


e(t)  = 


ue(t) 

i 


in  (15).  Balancing  the  forces  a  A  for  the  rod  with  those  of  the  restoring  mechanism  yields  the  lumped 
model 


Aod2uri+\  ,  CAdut  YA  d2U£  duri+\  1 

+  —rit]  +  —Ul(t)  =  ~  c‘n(t)  ~ kuc[t) 


+Aai[P(E(t))  -  Pr]  +  Aa2[P(E(t ))  -  PR] 2 


or,  equivalently, 

+  +  ku ^  =  d\[P{E{t))  -  Pr]  +  a2[P(E(t))  -  PR]2  (21) 

where 

CA  YA 

m  =  pA£  Emu  ,  c  =  —  +  q  ,  k=  -y  +  h  ,  a\  =  ^loi  ,  a2  =  ,4a2  (22) 

and  the  initial  conditions  are  ue(0)  =  uq  and  ^(0)  =  u\.  The  polarization  P  is  specified  by  the 
model  (12)  or  discretized  model  (14). 
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The  model  can  also  be  written  as  the  first-order  system 


ug(t)  =  Aue(t)  +  V(E(t)) 
ue(  0)  =  u0 


(23) 


where  ue(t)  =  [ug(t) ,  ug(t)]T ,  ug.{ 0)  =  [uo,ui]T  and 


A  = 


0  1 

—k/m  — c/m 


V(E(t))  =  1  [or (P(E(t))  -  PR)  +  a2(P(E(t ))  -  Pr)2]  J 


3.2  Cylindrical  Shell  Model 

To  quantify  the  dynamics  of  the  cylindrical  stage  depicted  in  Figure  2(b),  we  construct  a  linear 
shell  model  with  nonlinear  inputs  quantified  by  the  2-D  constitutive  relation  (16).  We  focus  on  the 
actuator  employed  for  transverse  displacements  since  real-time  control  of  this  component  is  required 
to  maintain  constant  forces  between  the  sample  and  nricro-cantilever.  The  mass  of  the  shell  employed 
for  horizontal  translation  is  combined  with  the  mass  of  the  sample  to  provide  an  inertial  force  acting 
on  the  free  end  of  the  vertical  actuator. 

For  modeling  purposes,  we  assume  that  the  shell  has  length  £,  thickness  h,  and  radius  R.  The 
axial  direction  is  specified  along  the  x-axis  and  the  longitudinal,  circumferential  and  transverse 
displacements  are  respectively  denoted  by  u,  v  and  w  as  depicted  in  Figure  8.  The  density  is  des¬ 
ignated  by  p  and  the  region  occupied  by  the  reference  or  middle  surface  of  the  shell  is  specified  by 
r0  =  [0,  i]  x  [0, 27 r].  In  accordance  with  the  constitutive  relations  (16),  Y,  C  and  v  denote  the  Young’s 
modulus,  Kelvin- Voigt  damping  coefficient  and  Poisson  ratio  for  the  material.  We  point  out  that 
ex,  eg  and  ex0  in  (16)  denote  strains  at  points  throughout  the  shell  thickness  whereas  2-D  shell  models 
are  formulated  in  terms  of  strains  ex,  eg  and  exg  in  the  reference  surface  of  the  shell.  The  relationship 
between  the  two  is  established  through  the  assumption  that  displacements  are  linear  through  the 
shell  thickness  which  comprises  one  of  the  fundamental  tenets  of  linear  shell  theory  [2,  25]. 

We  consider  the  case  in  which  the  bottom  edge  of  the  shell  (x  =  0)  is  clamped  and  the  opposite 
end  (x  =  £)  is  acted  upon  only  by  the  inertial  force  associated  with  the  combined  mass  m  of  the  x-y 
actuator  and  the  sample. 

As  detailed  in  [2,  25],  force  and  moment  balancing  yield  the  Donnell-Mushtari  shell  equations 


d2u  dNx 

Rp"m  ~  Rih 


dNxg 

de 


d2v  dNg  vdNxg 

Rph»  -  ~m  -  R^r  =  0 


(24) 


P  h—  _  pd2 Mx  _  1  d2Mg  n  Mxg 

ph  dt2  R  dx2  R  de2  2dxd0+  6 

where  the  force  and  moment  resultants  are 

Yh  Ch  h 

Nx  =  - 2^ ex  +  vee)  +  - +  v ^o)  -  — —  Wi(P  ~  Pr)  +  a2{P  ~  Pr )2] 

Yh  Ch  h 

Ng  =  l_i/2  (eg  +  vex)  +  1  _  ^2  (eg  +  vex)  -  ^[a i(P  -  PR)  +  a2(P  -  Pr)2] 


(25) 


Yh  Ch  . 

x9  ~  2(1  +  v) 6x9  +  2(1  +  0  6x0 
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u 


Figure  8:  Orientation  of  the  shell  geometry  used  when  quantifying  the  longitudinal,  circumferential 
and  transverse  displacements  u,  v  and  w. 


and 


Yh 3  Ch 3 

Mr  =  — - —{nx  +  vko)  +  tw - 2t(«x  +  vko) 


Me  = 


12(1  -  v2) 

Yh 3 

12(1  -  v2) 


(ko  +  vkx)  + 


12(1  -  v2) 
Ch 3 

12(1  -  v2) 


(ke  +  vkx) 


Yh?  Ch 3 

MxS  =  24(1  +  v)  Kxd  +  24(1  +  v) Kx6' 

The  midsurface  strains  and  changes  in  curvature  are 


du 

1  dv  w 

dv 

l  du 

e x  =  5 

Ox 

Ge=Rdd  +  R  ’ 

&x6  TT 

OX 

+  RdO 

d2  w 

i  d2w 

2  d2w 

dx 2 

’  Ke~  R2  d62 

i  ^x0  = 

R  dxdO 

The  boundary  conditions  for  the  fixed-end  at  x  =  0  are  taken  to  be 


(26) 


(27) 


dw 

u  =  v  =  w  =  — —  = 
ox 


0 


whereas  the  conditions 


Nt  =  —  m 


d2u 

W2 


„  l  dMx0  n 
®X  +  R  86  ~ 


Nxe  + 


Mxp 

~R~ 


=  0 


Mx  =  0 


are  employed  at  x  =  t.  The  first  resultant  condition  incorporates  the  inertial  force  due  to  the  mass 
m  of  the  PZT  actuator  employed  for  x-y  translation  along  with  the  mass  of  the  sample. 

To  reduce  smoothness  requirements  for  approximation  and  eliminate  the  Dirac  behavior  of  ex¬ 
ternal  inputs  at  x  =  l.  we  also  consider  a  weak  formulation  of  the  model.  The  state  is  taken  to  be 
z  =  (u(-,  -),u(£,  •))  in  the  state  space 


X  =  L2(D)  x  L2(D)  x  L2(D)  x  L2(0,2tt) 


where 

n  =  [0,£]  x  [0,2? r] 
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denotes  the  shell  region.  The  space  of  test  functions  is  specified  as 


v  =  {$  =  (0,  0,  <h,v)  £X\  01  6  Hl(p),<t>2  G  03  €  0?(0} 

where  rj(9)  =  0(00)  and 

H'(n)  =  {0eH\n)\(j)(o19)  =  o} 

(28) 

H${n)  =  {<!>€  H\n)  I  0(0,  9)  =  0(0,  0  =  0}  . 

Through  either  variation  principles  —  e.g.,  see  [2]  —  or  integration  by  parts,  one  obtains  the 
weak  formulation  of  the  thin  shell  model, 

L{Rph^h+ N‘a^+RN’°8t}dw=o  (29) 

i  -  rm*t£  -  -  5"*- TP  +  -  “• 

which  must  be  satisfied  for  all  $  6  V.  The  resultants  are  given  by  (25)  and  (26)  with  midsurface 
strains  and  changes  in  curvature  designated  in  (27). 

Remark  4  It  is  noted  that  the  d^i  poling,  used  to  generate  vertical  motion  in  the  stage,  produces  no 
polarization  contributions  to  the  moments.  However,  transverse  displacements  w  in  the  shell  model 
are  generated  by  the  Ng  resultant  in  the  w  relation  and  hence  all  three  components  of  the  displacement 
are  coupled. 


4  Model  Well-Posedness 


4.1  Rod  Model 


To  provide  a  framework  which  facilitates  the  establishment  of  criteria  which  guarantee  the  existence 
of  a  unique  solution  to  the  distributed  rod  model  with  nonlinear  inputs,  we  consider  a  Hilbert  space 
formulation  of  the  weak  model  formulation  (20)  with  the  state  and  test  function  spaces 

X  =  L2( 00)  x  R 

V  ={$  =  (</>,¥>)  exits  H\ 0,  0,  0(0)  =  0, 00  =  <0 


and  inner  products 


(4>i,<f>2)x=  /  pAefifcdx  +  mnpnp2 
Jo 

H 

($i,  $20  =  /  Y A<i>\ (i)2dx  +  kppi <p2 

Jo 


(30) 


where  <0  =  (0,  <^0, 0  =  (02, ^2)  with  ipx  =  0(0,  <p2  =  0(0- 

It  is  observed  that  V  is  densely  and  continuously  embedded  in  X  with  |$|x  <  c|4>|y;  this  is 
expressed  by  V  X.  Moreover,  when  one  defines  conjugate  dual  spaces  X*  and  V*  —  e.g.,  V* 
denotes  the  linear  space  of  all  conjugate  linear  continuous  functionals  on  k  —  two  observations  are 
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important:  (i)  X*  can  be  identified  with  X  through  the  Riesz  map,  and  (ii)  X*  V*.  Hence  the 
two  spaces  comprise  what  is  termed  a  Gelfand  triple  V  X  =  X*  c— >  V*  with  pivot  space  X  and 
duality  pairing  (duality  product)  (■, -)y«  V •  The  latter  is  defined  as  the  extension  by  continuity  of 
the  inner  product  (’>  ')x  from  V  x  X  to  V*  x  X.  Hence  elements  v*  6  V*  have  the  representation 
v*(v)  =  iv*,v)v*,v 

We  now  define  the  stiffness  and  damping  sesquilinear  forms  (Jt  :  V  x  V  — ►  C,  i  =  1,2,  by 


cri($i,  $2)  —  ($l,$2)y 

[l 

<t2($i,$2)=  /  CAfi^dx  A  cmpnp2- 


JO 

It  can  be  directly  verified  that  the  stiffness  form  satisfies 

(HI)  |ori($i,  $2)|  <  ci|<hi|y|$2|y  j  for  some  ci  G  M 

(H2)  ReiTi(<I>i,  <hi)  >  C2 1 d*  1 1 y-  ,  for  some  C2  >  0 

(H3)  ar($i,$2)  =  <h($2,$i) 

for  all  ip,<j)  EV.  Moreover,  the  damping  term  a2  satisfies 

(H4)  | cr2 (^1 ,  $2)|  <  C3 1 «J>i | v- 1 H>2 1 \x  ,  for  some  C3  G  M 

(H5)  Re  (72(^1,  $1)  >  C4 1 4*  1 1 y-  ,  for  some  C4  >  0 

The  input  space  is  taken  to  be  the  Hilbert  space  U  =  M  and  the 
defined  by 


(Bounded) 

(V-Elliptic) 

(Symmetric) 

(Bounded) 

(F-Elliptic). 
input  operator  B  :  U 


(31) 


(32) 
V*  is 


([H(E)](i),$)y,  y  =  [ai(P(E(t))  -PR)  +  a2(P(E(t))  -  PR )2]  [  <j>' dx  (33) 

Jo 

for  $  =  (1 j),  <p)  with  <p  =  4>(£).  It  is  observed  that  B  can  be  expressed  as 


[B(E)}(t)  =  [b(E)}(t)-g  ,geV* 

where 

[b(E)}(t)  =  (. P(E(t ))  -  PR)  +  a2(P(E(t ))  -  PR )2 

fi 

g{$)  =  <t>’  dx. 

Jo 

The  model  (20)  can  then  be  written  in  the  abstract  variational  formulation 
(u(t),  4>)y*y  +  Cr2(ll(t) ,  $)  +  0X  (u(t) ,  4>)  =  ([B(E)](t),  <E)V*,V 
n(  0)  =  no  ,  n(0)  =  u\ 

for  all  $  e  V. 

Alternatively,  one  can  define  the  operators  A\  £  C(V,  V*),i  =  1,2,  by 


(34) 


(35) 


(36) 


(Ai$  i,$2)y.,y  =  0*($1,$2) 
and  formulate  the  model  in  operator  form  as 

u(t)  +  A2u{t)  +  A\u{t)  =  [B(E)](t) 
n(  0)  =  no  ,  n(0)  =  ni 


(37) 


(38) 
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in  the  dual  space  V*.  This  formulation  illustrates  the  analogy  between  the  infinite-dimensional, 
strongly  damped  elastic  model  and  the  familiar  finite-dimensional  relations  (23). 

Model  Well-Posedness 

As  a  prelude  to  establishing  the  well-posedness  of  the  beam  model  with  hysteretic  E-P  relations, 
we  provide  a  lemma  which  quantifies  the  smoothness  of  the  input  operator. 

Lemma  1  Consider  field  inputs  E  £  (7[0,T].  The  input  operator  B  defined  by  (33)  then  satisfies 

B(E)  £  L2(0,  T;  V*).  (39) 

Proof.  In  Appendix  A,  we  establish  that  for  continuous  input  fields  E,  the  polarization  satisfies 
P  £  C[0,  T]  which  implies  that  b  defined  by  (35)  satisfies  &(•)  :  (7[0,T]  — >  C[0,T].  Hence  the  norm 

| [b(E)\ (t).  g(v)\ 

V *  =  sup  - 77 7j - 

vev  IMIv 

exists  for  each  t  £  [0,T],  Since  ||[R(.E,)](t)||vr*  =  |[&(F)](f)|  •  ||g||i/*,  it  follows  that 

l^(^)lli2(0,T;V)  <  {|[6(iT)](t)|2}  ■  T  ■  \\g\\v* 

which  implies  that 

B{E)  £  L2(0,  T;  V*). 

□ 

The  well-posedness  of  the  model  is  established  by  the  following  theorem  whose  proof  follows 
directly  from  Theorem  4.1  of  [2]  or  Theorem  2.1  and  Remark  2.1  of  [1]. 

Theorem  2  Let  a \  and  02  be  given  by  (31)  and  consider  continuous  field  inputs  E  £  (7[0,T].  There 
then  exists  a  unique  solution  w  to  (36),  or  equivalently  (38),  which  satisfies 

u  £  (7(0,  T;  V) 
u  £  (7(0,  T;  X). 

4.2  Shell  Model 

Similar  well-posedness  results  can  be  obtained  for  the  shell  model  (29)  through  consideration  of  an 
analogous  Hilbert  space  formulation  of  the  model.  Details  regarding  the  construction  of  appropriate 
inner  product  spaces,  sesquilinear  forms,  and  operators  can  be  found  in  [25,  30]. 

5  Numerical  Approximation  Techniques 

To  implement  the  distributed  models  for  either  the  rectangular  stacked  actuator  or  the  cylindrical 
actuator,  it  is  necessary  to  develop  appropriate  approximation  techniques  to  discretize  the  modeling 
PDE.  To  accomplish  this,  we  consider  general  Galerkin  methods  in  which  basis  functions  are  com¬ 
prised  of  spline  or  spline-Fourier  tensor  products.  The  resulting  methods  can  accommodate  a  variety 
of  boundary  conditions,  are  sufficiently  accurate  to  resolve  fine-scale  dynamics,  and  can  be  employed 
for  constructing  reduced-order  POD  approximates  for  real-time  implementation. 
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5.1  Stacked  Actuator  Model 

To  approximate  the  weak  form  of  the  stacked  actuator  model  (20),  we  employ  a  finite  element 
discretization  and  a  finite  difference  discretization  in  time.  The  semidiscrete  system  resulting  from 
the  finite  element  approximation  is  appropriate  for  finite  dimensional,  continuous  time  control  design 
whereas  the  fully  discrete  system  is  amenable  to  simulations  and  control  implementation. 

To  obtain  a  semidiscrete  system,  we  consider  a  uniform  partition  of  [0,£]  with  points  Xj  =  jh,j  = 
0, 1,  •  •  •  ,  N  with  step  size  h  =  £/N  where  N  denotes  the  number  of  subintervals.  The  spatial  basis 
{4>j}jLi  is  then  comprised  of  linear  splines 


( X  —  Xj- 1)  ,  Xj- 1  <  X  <  Xj 
(Xj. |_1  —  x)  ,  Xj  <  X  <  Xj. |_1 
0  ,  otherwise 


*  =  !,•••  ,  N  —  1 


1  j  (x  -  XAr_i)  ,  XN-!  <  x  <  xN 
MX)  =  h  0  .  otherwise 


(see  [21]  for  details  regarding  the  convergence  analysis  for  the  method).  The  solution  u(t,x )  to  (20) 
is  subsequently  approximated  by  the  expansion 

N 

UN{t,X )  =  ^Uj(t)(j)j{x)  . 

3= 1 

Through  construction,  the  approximate  solution  satisfies  the  essential  boundary  condition  uN(t,  0)  = 
0  and  can  attain  arbitrary  displacements  at  x  =  £. 

The  projection  of  the  problem  (20)  onto  the  finite  dimensional  subspace  VN  yields  the  semidis¬ 
crete  system 

z(t)  =  Az (t)  +  A  [ ai(P(t )  -  PR)  +  a2(P(t )  -  PR )2]  B 
z(0)  =  z0 

where  z (t)  =  [«i(£),  •  •  •  ,  URf(t),  «i(i),  •  •  •  ,  UAr(t)]T  and 


A  = 


0 

/tt —  In 


I 

/TT-lrf 


B  = 


0 

rxb 


The  mass,  stiffness  and  damping  matrices  have  the  components 


[M  }ij  =  { 


j  pA(f>i(f>j  dx  ,  i  /  N  or  j  /  N 

Jo 

re 

/  pA<f>i(j)j  dx  +  mi  ,  i  =  N  and  j  =  N 
-  Jo 


PK]ii  =  < 


o 

(  P 
1 0 


YAtfi'rf'j  dx  ,  N  or  j^N 


YAfiiftj  dx  +  ki  ,  i  =  N  and  j  =  N 


(40) 


(41) 
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and 


r 

/  cA^tp'j  dx  ,  i  /  N  or  j  /  N 

fe 

/  cA^cj)',  dx  +  C£  ,  i  =  N  and  j  =  N . 

.  io 


and  the  force  vector  is  defined  by 


[b]i  =  f  <t>'idx. 
Jo 


The  system  (40)  can  be  employed  for  finite-dimensional  control  design.  For  subsequent  implemen¬ 
tation,  we  consider  a  temporal  discretization  of  (40)  using  a  modified  trapezoid  rule.  For  temporal 
stepsizes  At,  this  yields  the  difference  equation 


zk+1  =  Wzfc  +  ^  aiP(tk)  +  aiP(tk+i)  +  a2P2(tk )  +  a2P2(4+i) 


VB 


(42) 


where  P  =  P  —  Pr,  tj  =  j'At,  z j  approximates  z (tj),  and 

-l 


W  =  (  I  -  —A 
2 


T  At 

H - A 

2 


¥  =  At  (  I  -  ^A 


-l 


This  yields  an  A-stable  method  requiring  moderate  storage  and  providing  moderate  accuracy. 


5.2  Cylindrical  Actuator  Model 

Due  to  the  inherent  coupling  between  longitudinal,  circumferential,  and  transverse  displacements  in 
combination  with  the  2-D  support  of  the  middle  surface,  the  numerical  approximation  of  the  model 
for  the  cylindrical  actuator  is  significantly  more  complicated  than  the  approximation  of  the  stacked 
actuator  model.  Among  the  issues  which  must  be  addressed  when  constructing  finite  element  or 
general  Galerkin  methods  for  the  shell  is  the  choice  of  elements  which  avoid  shear  and  membrane 
locking  and  the  maintenance  of  boundary  conditions.  We  summarize  here  a  spline-based  Galerkin 
method  developed  in  [7]  for  thin  shells  and  direct  the  reader  to  that  source  for  details  regarding  the 
construction  of  constituent  matrices  and  convergence  properties  of  the  method.  Details  regarding 
the  use  of  this  approximation  method  for  LQR  control  of  shells  utilizing  PZT  actuators  can  be  found 
in  [8]. 

The  bases  for  the  u,  v  and  w  displacements  are  respectively  taken  to  be 

$ufc(0,*)  =  zime  (t>Un{x)  ,  $Vk(d,x)  =  eime  (j)Vn{x)  ,  4v(<9,x)  =  eimd  </>Wn  (x) 

where  <!>Un,4>Vri  and  (f>Wn  are  cubic  B-splines  modified  to  satisfy  the  boundary  conditions  (e.g.,  see 
page  79  of  [21]).  The  approximating  subspaces  are 

f/f  =  span  {$uk}k=i  ,  VvN  =  span{$„j£1  ,  V*  =  span 

and  the  approximate  displacements  are  represented  by  the  expansions 

Nu 

uN(t,d,x)  =  y^uk(t)$Uk(0,x) 
k=  1 
Nv 

vN(t,  6,  x)  =  vk(t)$Vk  (9,  x )  (43) 

k= 1 
Nw 

wN(t,6,x)  =  y 'wk(t)$Wk(0,x). 
k= 1 
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The  restriction  of  the  problem  (29)  to  the  approximating  subspaces  and  construction  of  the 
forcing  vectors  subsequently  yields  the  matrix  system 

zN(t)  =  Az (t)  +  [ai(P(t)  -  Pr)  +  a2(P(t)  -  Pr)2}  B  ^ 

z(0)  =  z0  , 

where  z  =  #(i)]T,  with  t?(f)  =  [u(t),  v(i),  w(t)]T,  and 

0 

M_1b 


The  reader  is  referred  to  [7,  25]  for  details  concerning  the  construction  of  the  mass,  stiffness  and 
damping  matrices  M,  K  and  Q. 


6  Model  Validation 

Characterization  of  the  Stacked  Actuator 

We  consider  the  capability  of  the  modeling  framework  to  characterize  the  dynamics  of  the  stacked 
actuator  depicted  in  Figure  2(a).  The  PZT  actuator  had  a  length  of  £  =  2  x  10-2  nr  and  a  square 
cross-sectional  face  of  width  w  =  5  x  10-3  nr  so  that  the  cross-sectional  area  is  A  =  2.5  x  ICC5  nr2.  As 
illustrated  in  Figure  7,  one  end  of  the  actuator  was  considered  fixed  whereas  the  other  encountered 
elastic,  damping  and  inertial  effects  due  to  the  attached  components  of  the  stage  mechanism. 

To  validate  and  illustrate  properties  of  the  models,  we  consider  three  regimes:  (i)  end  displace¬ 
ments  quantified  by  the  lumped  model  (21)  with  the  thermally  inactive  kernel  (7)  employed  in  the 
polarization  model  (15),  (ii)  displacements  characterized  by  the  lumped  model  with  the  thermally 
active  polarization  kernel  (9),  and  (iii)  end  displacements  quantified  by  the  discretization  (42)  of  the 
distributed  model  (20).  It  is  illustrated  that  whereas  the  latter  choice  incorporates  the  distributed 
rod  nature  of  the  device,  the  fact  that  fields  and  stresses  are  uniform  along  the  rod  length  implies 
that  relative  displacements  are  also  uniform.  A  comparison  of  the  ODE  and  PDE  model  predictions 
at  the  rod  tip  (x  =  £)  illustrates  that  as  a  result,  the  ODE  provides  a  highly  accurate  characteriza¬ 
tion  with  significantly  less  computation  cost.  Hence  the  ODE  model  is  advantageous  for  real-time 
experimental  implementation. 

The  construction  of  the  models  requires  the  estimation  of  elastic,  damping  and  electromechanical 
parameters  in  addition  to  identification  of  the  densities  zq  and  v2  •  The  densities  were  estimated 
through  least  squares  fits  to  the  data  using  the  techniques  detailed  in  [25,  27].  The  manufacturer 
specifications  p  =  7600  kg/nr3  and  Y  =  7  x  1010  N/m2  were  employed  for  the  density  and  Young’s 
modulus  and  remaining  parameters  were  estimated  through  a  least  squares  fit  to  the  data.  The 
resulting  values  are  summarized  in  Table  1.  The  relation  between  the  rod  and  spring  parameters  is 
provided  by  (22). 

Lumped  Model  —  No  Thermal  Activation  in  Polarization  Relation 

We  consider  first  the  characterization  of  the  biased  minor  loop  data  shown  in  Figure  3  and 
frequency-dependent  data  from  Figure  4  using  the  lumped  model  (21)  with  the  thermally  inactive 
kernel  (7)  employed  in  the  polarization  model  (15).  It  should  be  noted  that  the  stage  was  disassembled 
between  the  quasistatic,  biased  minor  loop  experiments  and  the  frequency-dependent  experiments 
which  necessitated  the  re-identification  of  densities  for  the  two  cases. 
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Distributed  Model 

Parameter 

Value 

p  Y  C  mi  kn  ci  aq  <12 

7600  7  x  1010  5  x  106  4.015  8.49  x  10“5  440  1.54  x  1011  0 

Lumped  Model 

Parameter 

Value 

m  k  c  a\  ci2 

4.21  8.75  x  107  1.52  x  105  8.75  x  107  0 

Table  1:  Parameters  employed  in  the  distributed  (PDE)  model  (20)  and  lumped  (ODE)  model  (21) 
for  the  stacked  actuator. 


In  the  first  set  of  experiments,  displacement  data  measured  with  an  LVDT  was  collected  at  a 
sample  rate  of  0.1  Hz  and  four  input  field  levels  to  generate  a  set  of  biased  and  nested  transducer 
responses  ranging  from  nearly  linear  to  hysteretic  and  nonlinear  as  shown  in  Figures  3  and  9.  The 
densities  u\  and  v-i  and  parameters  summarized  in  Table  1  were  obtained  through  a  least  squares  fit 
to  the  full  data  set  comprised  of  four  loops.  The  resulting  model  accurately  quantifies  both  the  nest 
behavior  and  the  hysteresis  measured  at  increasing  input  levels. 

In  a  second  set  of  experiments,  data  was  collected  a  frequencies  ranging  from  0.279  Hz  to  27.9  Hz 
yielding  the  behavior  shown  in  Figure  4.  The  data  from  four  frequencies  was  used  to  re-identify 
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Figure  9:  Characterization  of  AFM  field-displacement  behavior  at  0.1  Hz  using  the  ODE  model  (21) 
with  the  thermally  inactive  kernel  (7). 
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parameters  in  the  polarization  model,  since  the  stage  had  been  modified,  thus  yielding  the  fits  shown 
in  Figure  10.  It  is  observed  that  the  model  characterizes  the  augmented  hysteresis  arising  at  higher 
frequencies  but  slightly  overpredicts  the  increase  in  displacement  following  field  reversal  which  is  due 
primarily  to  inertial  effects. 

Lumped  Model  —  Thermal  Activation  in  Polarization  Relation 

We  next  employ  the  thermally  active  kernel  (9)  in  the  polarization  model  to  incorporate  relaxation 
effects.  Parameters  in  the  polarization  model  were  again  identified  through  a  least  squares  fit  to  the 
four  frequency  data  sets  thus  yielding  the  model  fit  shown  in  Figure  11.  It  is  observed  that  use  of 
this  more  general  kernel  provides  additional  accuracy  at  higher  frequencies.  Whereas  this  improves 
characterization  capabilities,  the  added  accuracy  comes  at  the  cost  of  decreased  efficiency,  and  the 
criteria  of  accuracy  versus  efficiency  must  be  balanced  when  employing  the  model  for  real-time  control 
design  as  discussed  in  [12]. 


0  2000  4000  6000  8000 

Electric  Field  (V/m) 


0  2000  4000  6000  8000 

Electric  Field  (V/m) 


(a) 


(b) 


(c)  (d) 

Figure  10:  Characterization  of  AFM  field-displacement  behavior  using  the  ODE  model  (21)  with 
the  thermally  inactive  kernel  (7)  with  sample  rates  of  (a)  0.279  Hz,  (b)  1.12  Hz,  c)  5.58  Hz,  and 
(d)  27.9  Hz. 
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Figure  11:  Characterization  of  AFM  field-displacement  behavior  using  the  ODE  model  (21)  with  the 
thermally  active  kernel  (9)  with  sample  rates  of  (a)  0.279  Hz,  (b)  1.12  Hz,  c)  5.58  Hz,  and  (d)  27.9  Hz. 


Lumped  Versus  Distributed  Models 

It  has  been  observed  that  whereas  quantification  of  the  physics  of  the  stacked  actuator  leads  to 
the  rod  model  (20),  the  fact  that  stresses  and  fields  are  uniform  along  the  rod  length  implies  that 
relative  displacements  will  also  be  uniform.  This  motivates  consideration  of  the  lumped  model  (21) 
which  yielded  the  fits  shown  in  Figures  10  and  11. 

To  illustrate  the  validity  of  this  assumption,  the  difference  between  the  displacement  u(t,£), 
given  by  the  discretization  (42)  of  (20),  and  the  displacement  ug(t)  resulting  from  (21)  is  plotted  in 
Figure  12.  We  emphasize  that  when  constructing  the  PDE  model,  we  employed  the  parameter  values 
summarized  in  Table  1  which  are  consistent  with  the  spring  parameters  due  to  the  relation  (22).  The 
maximal  difference  of  5  x  10~10  is  5  orders  of  magnitude  less  than  the  micron-level  displacements 
being  characterized  thus  verifying  the  validity  of  the  ODE  model  in  this  regime.  The  accuracy  of 
the  ODE  model  has  important  ramifications  for  control  design  since  the  discretized  ODE  model  is 
significantly  more  efficient  to  implement  than  the  discretized  PDE  model. 
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Figure  12:  Difference  between  the  displacement  u(t,£)  given  by  the  distributed  model  (20)  and  Ui(t) 
given  by  the  lumped  model  (21). 


Characterization  of  the  Stacked  Actuator 

We  discuss  here  the  performance  of  the  cylindrical  shell  model  detailed  in  Section  3.2,  when  dis¬ 
cretized  using  the  Galerkin  techniques  summarized  in  Section  5.2,  for  characterizing  the  longitudinal 
displacements  of  the  cylindrical  PZT  shell  transducer  depicted  in  Figures  2(b)  and  7.  Whereas  the 
cylindrical  PZT  elements  employed  in  this  design  are  more  complex  than  the  rod  elements  used  in 
the  stage  design  depicted  in  Figure  2(a),  the  overall  transducer  is  simpler  and  has  the  advantage  of 
enhanced  vibration  isolation  and  diminished  hysteresis. 

The  experimental  cylindrical  transducer  had  a  length  of  £  =  0.0396  m,  radius  of  R  =  0.0056  nr 
and  thickness  h  =  0.0015  nr.  The  manufacturer  specifications  p  =  7600  kg/nr3  and  Y  =  7.1  x  1010 
N/nr2  were  employed  for  the  density  and  Young’s  modulus,  and  remaining  model  parameters  were 
estimated  through  a  least  squares  fit  to  the  data. 

The  longitudinal  displacement  uN  provided  by  (43)  is  compared  with  experimental  data  in  Fig- 
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Figure  13:  Characterization  of  the  relation  between  the  field  and  longitudinal  displacements  for  the 
cylindrical  actuator  depicted  in  Figures  2(b)  and  7. 
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ure  13.  We  note  that  due  to  the  inherent  coupling  between  the  longitudinal,  circumferential  and 
transverse  displacements,  u,v  and  w  in  the  model  (24),  the  approximate  displacements  (43)  are 
also  coupled  and  all  are  obtained  through  solution  of  (44)  —  we  plot  only  uN  since  it  corresponds 
to  measured  data.  The  nearly  linear  behavior  of  both  the  data  and  model  response  illustrates  an 
advantage  of  this  design  and  the  property  that  the  hysteretic  E-P  model  (12)  yields  approximately 
linear  behavior  in  low  drive  regimes.  The  fidelity  of  the  model  further  illustrates  the  accuracy  and 
flexibility  of  the  modeling  framework. 

7  Concluding  Remarks 

The  characterization  framework  developed  here  quantifies  both  the  approximately  linear  and  hys¬ 
teretic  properties  of  PZT  device  employed  in  atomic  force  microscope  (AFM)  positioning  mechanisms. 
In  the  first  step  of  the  development,  constitutive  relations  are  constructed  through  a  combination  of 
energy  analysis  at  the  lattice  level  and  stochastic  homogenization  techniques  based  on  the  assumption 
that  certain  parameters  are  manifestations  of  underlying  distributions.  These  relations  quantify  the 
frequency-dependent  hysteresis  exhibited  by  the  materials  for  general  drive  regimes  while  reducing 
to  approximately  linear  behavior  at  low  drive  regimes.  In  the  second  step  of  the  of  the  development, 
these  constitutive  relations  are  used  to  construct  lumped  and  distributed  rod  and  shell  models  for  the 
various  PZT  transducer  geometries.  The  accuracy  of  the  models  is  illustrated  through  comparison 
with  experimental  data  from  AFM  stages. 

An  important  property  of  the  framework  is  the  fact  that  resulting  models  can  be  approximately 
inverted  with  nearly  the  same  efficiency  as  the  forward  models  [12] .  This  provides  the  framework  with 
the  capability  for  providing  inverse  compensators  for  linear  control  design  [18,  19].  The  implemen¬ 
tation  of  feedback  control  designs  for  high  speed  scanning,  using  these  model-based  compensators, 
is  under  present  investigation. 

A  Continuity  of  the  Polarization  Model 

We  establish  here  the  continuity  of  the  homogenized  energy  model  (12), 

roo  roc 

[P(E)\(t)=  /  /  [P(E  +  Er,Ec,0}(t)u1(Ec)R2(EI)dEIdEc  (45) 

J  0  J — oo 

as  a  function  of  both  field  and  time  in  the  case  of  negligible  thermal  activation.  The  densities  v\  and 
i>2  satisfy  the  conditions  (13)  and  the  kernel  P  has  the  form 

P(E)  =  -  +  Pr6(E;Ec,Ei) 

V 

specified  in  (7). 

We  first  note  that  there  are  at  most  three  values  at  which  5  can  change  sign:  —Ec,  Ec  and 
— Ec  <  Et  <  Ec.  The  third  is  determined  by  the  initial  dipole  distribution  £,  as  depicted  in 
Figure  14(a),  and  is  typically  chosen  so  that  ER  =  0  when  E  +  Ej  =  0. 
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Figure  14:  (a)  Points  E  +  Ej  =  —Ec.  Ec  and  Et  at  which  5  =  ±1  changes  sign,  and  (b)  behavior  of 
5  associated  with  the  initial  dipole  distribution  at  E  +  Ej  =  Et- 


We  also  note  that  the  decay  conditions  (13)  dictate  that  v\  and  u2  satisfy  the  relations 


V2(Ei)\  <  C2 


i: 


u2(Ej)dE]  <  b2 


v\(Ec)dEc  <  bi 

where  b\,b2  and  c2  are  finite  constants. 

To  establish  the  continuity  of  P  with  respect  to  E,  we  consider  the  behavior  at  field  values  Eq 
and  Ei  where,  without  loss  of  generality,  we  take  Eq  <  E\.  When  integrating  with  respect  to  Ej ,  we 
decompose  the  interval  (—00,00)  into  seven  regions  delineated  by  the  points  —EC,EC,  Et  as  shown 
in  Figure  14(a).  For  this  decomposition,  we  note  that 

Eq)  ,  region  excludes 
—Ec,  Ec,  Et 

Eq)  +  2 Pr  ,  region  includes 
— Ec,  Ec,  Et  ■ 


f  1 


\P(Ei  +  Er,Ec,0-P(E0  +  Ei;Ec,0\  =  { 


To  consolidate  notation,  we  define  the  integrals 


I(a,b) 


-( E\  -  E0)u2(Ej)  dEj 
V 


IpR{a,b) 


(Ei  -  Eq)  +  2 Pr 


u2(Ej)  dEj. 
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It  subsequently  follows  that 


coo 

\P(Ei)  -  P(E0)\  <  /  {|/(-oo,  -Ec  -  Ei)\  +  \IPr{-Ec  -  Ei,-Ec  -  E0)\ 

Jo 

+1 1(—Ec  —  E0,  Et  —  Ei)\  +  | Ipr(Et  —  Ei,Et  —  E0 ) 
+| I{Et  ~  E0,EC  —  Ei)\  +  \IPr(Ec  -  Ei,Ec  -  £b)| 
+\I{EC  -  E0,oo)\}ui(Ec)dEc 


For  e  >  0,  take 


<  (Ei  —  Eq 


4c2 

- h  362 

V 


-(Ei  -  E0)  +  2 Pr 
V 


vi  (Ec)  dEc 


<  (Ei-E0)bi(—  +  3b2 


5  =  min 


V  V 


bi(*f  +  362 


-(Ei  -  Eq)  +  2 Pr 
V 


b(Ei  -  E0 )  +  2 PR 


Under  the  assumption  that  E  is  continuous  in  time  and  Eq  =  E(to),Ei  =  E(t\),  for  every  5  >  0 
there  exists  <5  >  0  such  that  if  \ti  —  to\  <  6,  we  are  guaranteed  that  \Ei  —  Eq\  <  5.  It  follows  that  if 
|6i  —  to  |  <  6)  the  polarization  values  satisfy  the  bound 


\[P(E)](ti)-[P(E)}(to)\<e 


thus  establishing  the  continuity  of  the  hysteresis  model.  This  holds  for  all  major  and  minor  loops.  As 
illustrated  in  Figure  6,  the  behavior  of  the  model  which  incorporates  thermal  activation  is  smoother 
than  the  thermally  inactive  case  considered  here.  For  brevity,  we  omit  the  proof  of  this  second  case. 
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